Kernel polynomial representation of imaginary-time Green's functions 
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Inspired by the recent proposed Legendre orthogonal polynomial representation of imaginary-time 
Green's functions, we develop an alternate representation for the Green's functions of quantum im- 
purity models and combine it with the hybridization expansion continuous-time quantum Monte 
Carlo impurity solver. This representation is based on the kernel polynomial method, which intro- 
duces various integral kernels to filter fluctuations caused by the explicit truncations of polynomial 
expansion series and improve the computational precision significantly. As an illustration of the new 
representation, we reexamine the imaginary-time Green's functions of single-band Hubbard model 
in the framework of dynamical mean-field theory. The calculated results suggest that with carefully 
chosen integral kernels the Gibbs oscillations found in previous orthogonal polynomial representa- 
tion have been suppressed vastly and remarkable corrections to the measured Green's functions have 
been obtained. 

PACS numbers: 71.10.Fd, 71.27.+a 



I. INTRODUCTION 

The rapid development of efficient numerical and an- 
alytical methods for solving quantum impurity models 
has been driven in recent years by the great success of dy- 
namical mean- field theory (DMFT) [l|-|3| and its non-local 
extensions. In the framework of DMFT, the momen- 
tum dependence of self-energy is neglected, then the so- 
lution of general lattice model may be obtained from the 
solution of an appropriately defined quantum impurity 
model plus a self-consistency condition. To solve the 
quantum impurity models numerous quantum impurity 
solvers have been developed in the last decades. [J Q In 
particular, continuous-time quantum Monte Carlo impu- 
rity solvers [7rfllj have become a very important tool for 
studying quantum impurity models, due to their accu- 
racy, efficiency and ability to treat extreme low tempera- 
ture and arbitrary interaction terms (for a recent review, 
see RefQjT 



Among various continuous-time quantum Monte Carlo 
algorithms, the hybridization expansion version (abbrevi- 
ated CT-HYB) 7-9] is the most powerful and reliable im- 
purity solver up to now and is widely used. In practice, as 
for the CT-HYB quantum impurity solver several severe 
technical limits still remain. One well-known problem is 
the high frequency noises commonly observed in the Mat- 
subara Green's function and the self-energy, [l| [l4| Sim- 
ilar problems also arise in the calculations of imaginary- 
time Green's functions and vertex functions, which are 
less emphasized in the literatures. In order to cure these 
problems, intuitive idea is to run more statistics in the 
Monte Carlo simulations to suppress the data fluctua- 
tions as far as possible. This strategy should mitigate 
the problems, but it will not avert them and will deteri- 
orate the efficiency of CT-HYB quantum impurity solver 
rapidly. Recently, Lewin Boehnke et a/. [151] have sug- 



gested to measure the single-particle and two-particle 
Green's function in the basis of Legendre orthogonal 
polynomials. The orthogonal polynomial method (OPM) 
provides a more compact representation of Green's func- 
tions than standard Matsubara frequencies, and there- 
fore significantly reduces the memory storage size of these 
quantities. Moreover, it can be used as an efficient noise 
filter for various physical quantities within the CT-HYB 
quantum impurity solver: the statistical noise is mostly 
carried by high-order Legendre coefficients which should 
be truncated, while the physical properties are deter- 
mined by low-order coefficients which should be retained. 
By and by the OPM is used for the computations of 
single-particle Green's function and lattice susceptibil- 
ities in the context of realistic DMFT calculations in 
combination with the local density approximation to the 
density functional theory (LDA+DMFT) . [ll| [l?} By us- 
ing the Legendre orthogonal polynomial representation, 
the accuracy of CT-HYB impurity solver is greatly im- 
proved. But according to careful examinations, the so- 
called Gibbs oscillations can be easily found in the re- 
sulting Green's functions and other physical quantities, 
which may be mainly due to the rough truncation of 
Legendre basis. The sign of Gibbs oscillations is that 
the resulting physical quantities are smooth but oscillat- 
ing periodically with the scattered direct measurements. 
The situation is even worse in the insulating state where 
the Gibbs oscillations will cause the reconstructed single- 
particle Green's function to break the causality. Noted 
that a common procedure to damp these oscillations re- 
lies on an appropriate modification of the expansion coef- 
ficients by some integral kernels, which is the well-known 
kernel polynomial representation. [l8| Thus we adopt the 
kernel polynomial method (KPM) to improve the mea- 
surement of single-particle and two-particle quantities 
within CT-HYB quantum impurity solver and expect to 
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obtain significant improvements. 

The rest of this paper is organized as follows: In 
Sec lII Al a brief introduction to the orthogonal polynomial 
representation is provided. The original implementation 
is based on Legendre polynomials only, and a straight- 
forward generalization to Chebyshev polynomials is pro- 
posed. In Sec lII Bl the kernel polynomial representation 
is presented in details. Then in Sec lIIII we benchmark the 
KPM by reexamining the imaginary-time Green's func- 
tion of single-band half-filled Hubbard model, and the 
characteristics of different integral kernel functions which 
are used to alter the expansion coefficients are discussed. 
Section IIVI serves as a conclusion and outlook. Finally 
in Appendix [A] concise introductions for the Chebyshev 
and Legendre orthogonal polynomial series are available 
as well. 



T > 
T < 0. 



(6) 



Lewin Bochnke et al. [15|] have chosen the Legendre or- 
thogonal polynomials as their preferred basis to expand 
single-particle and two-particle Green's functions. But 
it should be stressed that a priori different orthogonal 
polynomial bases may be used as well. Thus we try to 
generalize the OPM to use Chebyshev orthogonal poly- 
nomials as an optional basis. It is well-known that there 
exist two kinds of Chebyshev polynomials. By using 
the Chebyshev polynomials of second kind U n (x) as ba- 
sis, the imaginary-time Green's functions G(r) can be 
expressed by the following equations, 



II. METHOD 



A. Orthogonal polynomial representation 



' n>0 



U n {x( T )]G r , 



(7) 



In the OPM, the imaginary-time Green's function G(r) 
where r G [0, ft] can be expanded in terms of Legen- 
dre orthogonal polynomials P n (x) defined on the interval 
[-1,1], 



i 

G ( T ) = 8 V2^+lP n [x(r)]G n , 



n>0 



G n = V^TT / drP„[x(r)]G(r) 
Jo 



(1) 



(2) 



where ft is inverse temperature, x(t) = =£■ — 1 and G n 
denotes the expansion coefficients of G(r) in the Legen- 
dre orthogonal polynomials basis. 15] Since the expansion 
coefficients generally show a very fast decay with n, the 
expansion in Legendre polynomials can be truncated at 
a maximum order n max . In the CT-HYB quantum im- 
purity solver, the formula for measuring imaginary-time 
Green's function G(r) @-||| is 



Ik k 



= 1 3 = 1 



A(r,r') 



5(t - t'), t' > 

-S(T-r'-ft), r'<0, 



(3) 



(4) 



where k is the order of diagrammatic perturbation expan- 
sion series, matrix element (M" 1 )^ = F{rf — rj) where 
F(t) is the hybridization function, rf and r| are the co- 
ordinates in imaginary-time axis for create and destroy 
operators, respectively. By utilizing Eq.Q and Eq.©, 
the Legendre coefficients for G(r) finally become 



G,. 



V2n + 1 

ft 



/ k k 
\i=l j=l 



(5) 



G n = - dTU n [ X {T)]y/l-x(T)*G(T). (8) 



After a straightforward substitute, the Chebyshev coeffi- 
cients for G(r) finally become 



k k 



where 



Un(r) 



Un[x(T)l 



T > 



and 



x{t) 



-U n [x{T + ft% T<0, 



x(t), t > 

x(r + ft), t < 0. 



(9) 
(10) 



(11) 



The CT-HYB quantum impurity solver can directly ac- 
cumulate the Legendre or Chebyshev coefficients G n in- 
stead of original Green's functions G(r). Once the Monte 
Carlo sampling has been finished, G(r) can be recon- 
structed analytically by using the expansion coefficients. 
Since the coefficients decay very quickly, the orthogo- 
nal polynomial bases are much more compact and are 
particularly interesting for storing and manipulating the 
two-particle quantities, like vertex function etc. Further- 
more, the Monte Carlo noises are mainly concentrated in 
the high-order expansion coefficients, and the numerical 
values of them are usually very small. So a rough trun- 
cation method can be developed to filter out the noises 
and obtain more smooth and accurate results. 
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TABLE I. Summary of different integral kernel functions f n that can be used to improve the quality of an order N Chebyshev 
or Legendre series 



name 


fn 


parameters positive 


Jackson 


Tf [{N - n + 1) cos^) + sin(^_) cot( 1 ^ T )_ 


none 


yes 


Lorentz 


sinh[A(l - n/JV)]/ sinh(A) 


A G n 


yes 


Fejer 


1-n/N 


none 


yes 


Wang-Zunj 


;er exp [- (a§)^ 


a, pen 


no 


Dirichlet 


1 


none 


no 
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FIG. 1. (Color online) Classic integral kernels f n used to 
improve the quality of polynomial expansion series. In this 
figure, the order of expansion series is N = 64. The Dirichlet, 
Jackson, and Fejer kernels take no parameters. For Lorentz 
kernel, the A parameter is fixed to be 1.0. And for Wang- 
Zunger kernel, the a and f3 parameters are 1.0 and 4.0, re- 
spectively. 



B. Kernel polynomial representation 



alent to previous OPM. In addition to the Dirichlet ker- 
nel, other classic integral kernel functions, like Jackson, 
Lorentz, Fejer, and Wang-Zunger etc., are collected in 
Tab U and plotted in Fig Irrespectively. Note that for all 
the kernels fo must be equal to 1 and /i must approach 
1 as n — > oo. The optimal integral kernel function par- 
tially depends on the considered application. According 
to the literature, [l8[ the Jackson kernel may be the best 
for most applications, the Lorentz kernel may be the best 
for Green's function, while the Fejer kernel is mainly of 
academic interest. 

Finally, we note that the integral kernel functions f n 
can be evaluated and stored in advance, so the KPM has 
not effect on the computational efficiency of CT-HYB 
quantum impurity solver. The implementation of KPM 
is very simple, only small modifications are needed for the 
OPM's version of CT-HYB, i.e., replacing G„ by G n f n . 
Since the kernel polynomial and orthogonal polynomial 
representations are only alternate bases for single-particle 
and two-particle quantities, so both methods can be im- 
plemented in segment picture. 7] and general matrix Q 
formulations of CT-HYB impurity solvers to improve the 
accuracy and efficiency. 



The basic idea of OPM is to expand single-particle 
Green's functions G(r) in infinite series of Chebyshev or 
Legendre polynomials, and then use Monte Carlo algo- 
rithm to sample the expansion coefficients G n directly. 
As expected for a numerical approach, however, the ex- 
pansion series will remain finite actually, and we thus 
arrive at a classical problem of approximation theory. In 
our case the problem is equivalent to find the best ap- 
proximation to G(r) given a finite number of G n . Expe- 
rience shows that a simple truncation of the infinite se- 
ries leads to poor precision and fluctuations, which also 
known as Gibbs oscillations. [l8| For examples, as for the 
reconstructed Green's function G(r) in insulating state, 
almost periodic Gibbs oscillations are clearly identified 
in a wide r range. 

A common procedure to damp the Gibbs oscillations is 
to introduce some kind of integral kernel function f n and 
change the expansion coefficients from G n to G„/ ra .[l8[ 
Obviously, the simplest integral kernel function, which 
is usually attributed to Dirichlet, is obtained by setting 
/„ = 1. By using the Dirichlet kernel, the KPM is equiv- 



III. BENCHMARK 

In this section, we try to benchmark the kernel poly- 
nomial representation and compare the calculated results 
with those obtained by orthogonal polynomial represen- 
tation and conventionally direct measurements. For the 
sake of simplicity, a single-band Hubbard model on Bcthc 
lattice is used as a toy model to examine our implemen- 
tations of OPM and KPM. Here U = 4.0 and /3 = 10.0 
for metallic case, and U — 6.0 and (3 — 50.0 for insulating 
case. The band is with bandwidth 2.0, and a semicircular 
density of states is chosen. The chemical potential fi is 
fixed to be U/2 to keep the model under half- filling. Un- 
less it is specifically stated, this model is used throughout 
this section. This toy model is studied in the framework 
of single site DMFT0, 0] and the segment picture ver- 
sion of CT-HYB is used as quantum impurity solver. 
In each DMFT iterations, typically 4 x 10 8 Monte Carlo 
sweeps have been performed to reach sufficient numerical 
precision. 
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FIG. 2. (Color online) Chebyshev and Legendre coefficients 
/3\G n \ of imaginary-time Green's functions of the single-band 
half-filled Hubbard model on the Bethe lattice within DMFT. 
The Coulomb interaction strength U is 4.0 and inversion tem- 
perature /3 is 10. The coefficients in the pink region contribute 
very little to the resulting Green's function. 



A. Metallic state 



Let's first concentrate our attentions to the metal- 
lic state. In figure [5] the "bare" expansion coefficients 
(3\G n \ of Chebyshev and Legendre orthogonal polyno- 
mials are shown. As is pointed out by Lewin Boehnke 
et aZ. [ToT| . due to the constraint of particle- hole symme- 
try the expansion coefficients for odd order n should be 
zero. Indeed, the coefficients in our data for odd n's all 
take on a very small value, compatible with a vanish- 
ing value within their error bars. The even n coefficients 
instead show a very fast decay. As is shown in this fig- 
ure, n max = 15 ~ 25 is enough for both Chebyshev and 
Legendre polynomial representations to obtain converged 
and accurate results. In our simulations, n max is fixed to 
be 24. 

Figure [3] shows the calculated imaginary-time Green's 
function G(t) by using KPM with different orthogonal 
polynomials and integral kernel functions. It is appar- 
ent that the directly measured G(t) is full of noises and 
fluctuations, which are neg ative for the later analytical 
continuation procedure. |2fj| Once the OPM is used (i.e., 
Dirichlet kernel /„ = 1 is adopted), G(r) turns smooth 
but obvious undulations still exist. If the Lorentz, Fejer, 
and Wang-Zunger kernels are applied one by one, the 
Green's functions are smooth and without obvious un- 
dulations, but deviate systematically from the scattered 
data. As a general view, the Jackson kernel function 
is the optimal choice. The resulting G(r) evaluated by 
Jackson kernel function is smooth and nicely interpolates 
the directly measured data. As is expected, the type of 
orthogonal polynomials has little impact to the interpo- 
lated G(t). It seems that the Chebyshev polynomials do 
a bit better than Legendre polynomials. 
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FIG. 3. (Color online) The imaginary-time Green's function 
G(t) in t e [0.6/?, 0.8/3] interval of the single-band half-filled 
Hubbard model. The Coulomb interaction strength U is 4.0 
and /3 = 10. Upper panel: G(r) calculated by Chebyshev 
polynomials with or without integral kernel functions. Lower 
panel: G(t) calculated by Legendre polynomials with or with- 
out integral kernel functions. The order for polynomial expan- 
sion series is 24. 



B. Insulating state 

Now let's turn to the insulating state. When U = 6.0 
and /3 = 50 an definitely insulating solution is obtained 
within DMFT. The "bare" Chebyshev and Legendre co- 
efficients /3|G„| of G(t) are shown in Figj4] Just similar 
to the metallic state, G„ takes very small value for odd n 
and can be ignored safely, and for even n, G n converges to 
zero very quickly. For Chebyshev and Legendre polyno- 
mials, n max = 35 ~ 45 or n max = 40 ~ 50, respectively. 
Thus the Chebyshev polynomials is more compact and 
efficient than Legendre polynomials. In current simula- 
tions, n max is fixed to be 64 uniformly. 

The calculated imaginary-time Green's function G(r) 
by using KPM with different orthogonal polynomials and 
integral kernel functions are illustrated in FigJSJ A cur- 
sory look could lead you to believe that the reconstructed 
Green's functions by KPM or OPM agree rather well with 
the scattered data. Next let's zoom in r € [0.2/3, 0.8/3] 
interval and check carefully the magnified G(r), which 
are just depicted in the insets of FigJS] Clearly, G(r) 
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FIG. 4. (Color online) Chebyshev and Legendre coefficients 
/3\G n \ of imaginary-time Green's functions of the single-band 
half-filled Hubbard model on the Bethe lattice within DMFT. 
The Coulomb interaction strength U is 6.0 and inversion tem- 
perature ft is 50. The coefficients in the pink region contribute 
very little to the resulting Green's function. 
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is very close to zero in this region. The scattered data 
obtained by direct measurement exhibit periodical undu- 
lations. The reconstructed G(t) by OPM (with Dirichlet 
kernel) displays stronger periodical oscillations and vio- 
lates the causality at the same time, which means that 
the results will be even deteriorated by using orthogo- 
nal polynomial representation. The results obtained by 
Wang-Zunger kernel fit original data very well and obvi- 
ous improvement is not observed. As for the Lorentz and 
Fejer kernels, the calculated results deviate the scattered 
data systematically. Again, the Jackson kernel is the op- 
timal choice. The calculated results are very smooth, 
perfectly interpolate the scattered data, and obey the 
causality. 



IV. CONCLUSIONS 

It is suggested that the OPM based on Legendre or- 
thogonal polynomials can be used to improve the ac- 
curacy and computational efficiency of CT-HYB quan- 
tum impurity solver. [l5l| In this paper, we develop a 
better representation to calculate the single-particle and 
two-particle quantities. Firstly, we generalize the OPM 
to Chebyshev orthogonal polynomial basis. Secondly, 
the KPM based on various integral kernel functions is 
proposed to damp the Gibbs oscillations observed in 
the single-particle and two-particle Green's functions ob- 
tained with OPM and improve the accuracy of them fur- 
ther. According to the benchmark results for single-band 
half-filled Hubbard model, it is demonstrated that the 
Jackson kernel is the optimal choice for imaginary-time 
Green's function G(r) and other quantities. Though the 
KPM presented in this paper is mainly developed for the 
CT-HYB quantum impurity solver, it can be easily gen- 
eralized to other continuous-time quantum Monte Carlo 
impurity solvers. [12] 
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Appendix A: Chebyshev and Legendre polynomials 



FIG. 5. (Color online) The imaginary-time Green's func- 
tion G(r) of the single-band half-filled Hubbard model. The 
Coulomb interaction strength U is 6.0 and /3 = 50. Upper 
panel: G(t) calculated by Chebyshev polynomials with or 
without integral kernel functions. Lower panel: G(t) calcu- 
lated by Legendre polynomials with or without integral kernel 
functions. The insets in both panels show the fine structures 
of G(t) in r 6 [0.2/3, 0.8/3] interval. The order for polynomial 
expansion series is 64. 



In this appendix, we summarize for convenience some 
basic properties of the Chebyshev and Legendre polyno- 
mials, and the first few Chebyshev and Legendre polyno- 
mials are illustrated in Figj6l Further reference can be 
found in Refill 

The Chebyshev polynomials are a sequence of orthog- 
onal polynomials which can be defined recursively. One 
usually distinguishes between Chebyshev polynomials of 
the first kind which are denoted by T n and Chebyshev 
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FIG. 6. (Color online) Chebyshev and Legendre orthogonal 
polynomials. Upper panel: The first few Chebyshev polyno- 
mials of the second kind U n in the domain 1 < x < 1. Lower 
panel: The first few Legendre polynomials P n in the domain 
1 < x < 1. 



polynomials of the second kind which are denoted by U n . 
The Chebyshev polynomials T n or U n are polynomials of 
degree n. T n and U n are defined by the following recur- 
rence relations: 

T (x) = 1, (Al) 
Ti(x)=x, (A2) 
T n+ i{x) = 2xT n (x) - T n -i(x). (A3) 

and 

U (x) = 1, (A4) 



Ui(x) = 2x, 



(A5) 



U n+1 (x) = 2xU n (x) - U n -x(x). (A6) 
The Chebyshev polynomials of the first and second kind 
are closely related by the following equations: 



and 



T n+1 (x) = xT n {x) - (1 - ir 2 )C/„_i(x), (A7) 



T n (x) = U n (x) - x U n -!(x). (A8) 



The Chebyshev polynomials of the first kind are orthog- 



onal with respect to the weight 
[—1,1], i.e. we have: 



on the interval 



T^n (x)T m (x) 



dx 



VT - 



: n 7^ m, 
= ^ 7r : n = m = 0, (A9) 
7r/2 : n = m ^ 0. 



Similarly, the Chebyshev polynomials of the second kind 
are orthogonal with respect to the weight Vl — x 2 on the 
interval [1, 1], i.e. we have: 







U n (x)U m (x)Vl-X 2 dx = i" ■~r-, (A1Q) 

I 7T/2 : n — m. 



Similar to the Chebyshev polynomials, Legendre poly- 
nomials P n are orthogonal polynomials of degree n, which 
can be defined by the following recurrence relations: 



P Q (x) = 1, 



Pi (a;) = x, 



(All) 
(A12) 



(n + l)P„+i(a:) = (2n + l)xP n {x) - nP n -i(x). (A13) 

The orthogonality of Legendre polynomials on the inter- 
val [-1,1] reads: 



P m (x)P n (x) dx 



2n + 1 



(A14) 
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